{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ccd542f6",
   "metadata": {},
   "source": [
    "# 滤波器实现"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c0d1e46",
   "metadata": {},
   "source": [
    "## 卷积和滤波"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f82108fa",
   "metadata": {},
   "source": [
    "滤波的数学基础是卷积。对于有限脉冲响应(FIR)滤波器，滤波运算的输出$y(k)$是输入信号$x(k)$与脉冲响应$h(k)$的卷积：\n",
    "$$ y(k)=\\sum_{l=-\\infty}^{\\infty}h(l)x(k-l) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28deba60",
   "metadata": {},
   "source": [
    "如果输入信号也是有限长度的，您可以使用conv/convolve函数来执行滤波运算。例如，要用三阶平均值滤波器对包含五个样本的随机向量进行滤波，可以将$x(k)$存储在向量x中，将$h(k)$存储在向量h中，并求这两个向量的卷积："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0a10fef0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "\n",
    "x = np.random.randn(5)\n",
    "h = np.ones(4)/4\n",
    "y = signal.convolve(x,h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c478b200",
   "metadata": {},
   "source": [
    "y的长度比x和h的长度之和小1。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c3471cb",
   "metadata": {},
   "source": [
    "## 滤波器和传递函数"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04a86c69",
   "metadata": {},
   "source": [
    "滤波器的传递函数是其脉冲响应的Z变换。对于FIR滤波器，输出y的Z变换$Y(z)$是传递函数和输入x的Z变换$X(z)$的乘积："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "886ec3a1",
   "metadata": {},
   "source": [
    "$$Y(z)=H(z)X(z)=\\left(h(1)+h(2)z^{-1}+\\cdots+h(n+1)z^{-n}\\right)X(z)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b5dc22",
   "metadata": {},
   "source": [
    "多项式系数$h(1),h(2),...,h(n+1)$对应于第n阶滤波器的脉冲响应的系数。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07647b5a",
   "metadata": {},
   "source": [
    "**注意**：Matlab中滤波器系数索引从1到(n + 1)，而Python中是从0到n。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e8e1f2e",
   "metadata": {},
   "source": [
    "FIR 滤波器也称为全零、非递归或移动平均值(MA)滤波器。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1349b5a2",
   "metadata": {},
   "source": [
    "对于无限脉冲响应(IIR)滤波器，传递函数不是多项式，而是有理函数。输入和输出信号的Z变换的关系是："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d539bec",
   "metadata": {},
   "source": [
    "$$Y(z)=H(z)X(z)=\\frac{b(1)+b(2)z^{-1}+\\cdots+b(n+1)z^{-n}}{a(1)+a(2)z^{-1}+\\cdots+a(m+1)z^{-m}}X(z)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1b0d77c",
   "metadata": {},
   "source": [
    "其中$b(i)$和$a(i)$是滤波器系数。在本例中，滤波器的阶是n和m的最大值。n = 0的IIR滤波器也称为全极点、递归或自回归(AR)滤波器。n和m均大于零的IIR滤波器也称为极点-零点、递归或自回归移动平均值(ARMA)滤波器。缩写AR、MA和ARMA通常应用于与滤波随机过程相关联的滤波器。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d7f8952",
   "metadata": {},
   "source": [
    "## 使用filter函数进行滤波"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "391dfd47",
   "metadata": {},
   "source": [
    "对于IIR滤波器，滤波运算不能用简单的卷积来说明，而需要用可从传递函数关系中找到的差分方程来说明。假设$a(1)=1$，将分母移到左侧，并进行逆Z变换，以获得"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7e63b12",
   "metadata": {},
   "source": [
    "$$y(k)+a(2)y(k-1)+\\cdots+a(m+1)y(k-m)=b(1)x(k)+b(2)x(k-1)+\\cdots+b(n+1)x(k-n)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8406d875",
   "metadata": {},
   "source": [
    "对于当前输入、过去的输入以及过去的输出，$y(k)$是"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc9e60e7",
   "metadata": {},
   "source": [
    "$$y(k)=b(1)x(k)+b(2)x(k-1)+\\cdots+b(n+1)x(k-n)-a(2)y(k-1)-\\cdots-a(m+1)y(k-m)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f46c09b",
   "metadata": {},
   "source": [
    "这是数字滤波器的标准时域表示。从$y(1)$开始，假设一个初始条件为零的因果系统，表示形式等效于"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38d82159",
   "metadata": {},
   "source": [
    "$y(1)=b(1)x(1)$  \n",
    "$y(2)=b(1)x(2)+b(2)x(1)-a(2)y(1)$  \n",
    "$y(3)=b(1)x(3)+b(2)x(2)+b(3)x(1)-a(2)y(2)-a(3)y(1)$  \n",
    "$\\quad \\; \\; \\;\\vdots$  \n",
    "$y(n)=b(1)x(n)+\\cdots+b(n)x(1)-a(2)y(n-1)-\\cdots-a(n)y(1)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5ffe10e",
   "metadata": {},
   "source": [
    "要实现此滤波运算，您可以使用MATLAB的filter函数或Python的lfilter函数。filter/lfilter将系数存储在两个行向量中，其中一个为分子，一个为分母。例如，要求解差分方程"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dea31361",
   "metadata": {},
   "source": [
    "$$y(n)-0.9y(n-1)=x(n)\\quad \\Rightarrow \\quad Y(z)=\\frac{1}{1-0.9z^{-1}}X(z)=H(z)X(z)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fa86cf6c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "b = np.array([1])\n",
    "a = np.array([1,-0.9])\n",
    "zi = signal.lfilter_zi(b,a)*0\n",
    "y,_ = signal.lfilter(b,a,x,zi=zi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c31a8a4",
   "metadata": {},
   "source": [
    "filter/lfilter提供的输出样本的数量与输入样本一样多，即y的长度与x的长度相同。如果a的第一个元素不是1，则filter/lfilter在实现差分方程之前，将系数除以$a(1)$。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
